R jako GIS

Podpůrný text pro předmět: GIS pro biologické aplikace
Autor: Matěj Man
Aktualizace: 18. 11. 2019
Odkaz na R skript: Kopírovat a vložit do R

Další zdroje:

Knihovny

list.of.packages <- c("sf","raster","mapview","whitebox","randomcoloR","leaflet")
# install.packages(list.of.packages)

library(sf)
library(raster)
library(mapview)
library(randomcoloR)
library(leaflet)

Načítání vektorových GIS dat do R

## Načítání vektorových dat shp
# nastavte kde máte u sebe na PC data
# pozor musíte zdvojit nebo otočit lomítka
cesta<-"d:\\Owncloud\\ŠKOLA\\Učení\\GIS\\cv 6\\2020_2021_cv6\\" 
# zkonstuuje cestu kde leží vrstva Brdy
data.path<-paste0(cesta,"\\parky2\\CHKO_Brdy.shp")
# načte .shp do R
brdy<-st_read(data.path,stringsAsFactors = F,quiet = T)
# zobrazí načtený soubor v interkativním okně mapy
# mapview(brdy)

# obdobně načteme třeba hranici ČR
data.path<-paste0(cesta,"\\hrcr1_wgs.shp")
hrcr<-st_read(data.path,stringsAsFactors = F,quiet = T)

# prostý obrázek, bez interaktivity
plot(st_geometry(hrcr)) 
# parametr add přidá další vrstvu do existujícího obrázku
plot(st_geometry(brdy),col="red",add=T) 

# Nastaví pracovní adresář working dorectory (WD)
# dále nemusíme data z WD volat celou cestou, stačí názvy 
setwd(cesta) 
## Načítání vektorových dat z tabulky
# načíst prostou tabulku
chmu<-read.table("staniceCHMUtablecoma.csv",header = T, sep=",")
# převést tabulky na prostorová data
chmu_sf<-st_as_sf(chmu, coords = c("Xcoo", "Ycoo"),crs = 4326) 

## vizualizovat prostorová data interaktivně
# pokročilejší balíček leaflet

# definujeme paletu o 17 barvách podle typu stanice
distCol<- colorFactor(distinctColorPalette(17), chmu_sf$Station.ty)

# definoujeme okno s mapou 
leaflet(chmu_sf) %>% 
  addTiles() %>%  
  addCircleMarkers(popup=chmu_sf$Name,color = ~distCol(Station.ty),
                   stroke = FALSE, fillOpacity = 0.8,radius=4) %>%
  addLegend(pal = distCol, values = ~Station.ty)

Projekce a crs

st_crs(brdy) # jaký crs má vrstva brdy ?
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
st_crs(hrcr) # jaký crs má vrstva hrnice čr ?
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
## vektory
# transformace podle EPSG
brdy_32633<-st_transform(brdy,32633)
st_crs(brdy_32633) # jaký crs má vrstva brdy_32633 ?
## Coordinate Reference System:
##   User input: EPSG:32633 
##   wkt:
## PROJCRS["WGS 84 / UTM zone 33N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 33N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",15,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World - N hemisphere - 12°E to 18°E - by country"],
##         BBOX[0,12,84,18]],
##     ID["EPSG",32633]]
# transformace podle proj4 string
brdy_5514<-st_transform(brdy, "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56")
st_crs(brdy_5514) # jaký crs má vrstva brdy_5514 ?
## Coordinate Reference System:
##   User input: +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 
##   wkt:
## BOUNDCRS[
##     SOURCECRS[
##         PROJCRS["unknown",
##             BASEGEOGCRS["unknown",
##                 DATUM["Unknown based on Bessel 1841 ellipsoid",
##                     ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
##                         LENGTHUNIT["metre",1,
##                             ID["EPSG",9001]]]],
##                 PRIMEM["Greenwich",0,
##                     ANGLEUNIT["degree",0.0174532925199433,
##                         ID["EPSG",9122]]]],
##             CONVERSION["unknown",
##                 METHOD["Krovak (North Orientated)",
##                     ID["EPSG",1041]],
##                 PARAMETER["Latitude of projection centre",49.5,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8811]],
##                 PARAMETER["Longitude of origin",24.8333333333333,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8833]],
##                 PARAMETER["Co-latitude of cone axis",30.2881397222222,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",1036]],
##                 PARAMETER["Latitude of pseudo standard parallel",78.5,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8818]],
##                 PARAMETER["Scale factor on pseudo standard parallel",0.9999,
##                     SCALEUNIT["unity",1],
##                     ID["EPSG",8819]],
##                 PARAMETER["False easting",0,
##                     LENGTHUNIT["metre",1],
##                     ID["EPSG",8806]],
##                 PARAMETER["False northing",0,
##                     LENGTHUNIT["metre",1],
##                     ID["EPSG",8807]]],
##             CS[Cartesian,2],
##                 AXIS["(E)",east,
##                     ORDER[1],
##                     LENGTHUNIT["metre",1,
##                         ID["EPSG",9001]]],
##                 AXIS["(N)",north,
##                     ORDER[2],
##                     LENGTHUNIT["metre",1,
##                         ID["EPSG",9001]]]]],
##     TARGETCRS[
##         GEOGCRS["WGS 84",
##             DATUM["World Geodetic System 1984",
##                 ELLIPSOID["WGS 84",6378137,298.257223563,
##                     LENGTHUNIT["metre",1]]],
##             PRIMEM["Greenwich",0,
##                 ANGLEUNIT["degree",0.0174532925199433]],
##             CS[ellipsoidal,2],
##                 AXIS["latitude",north,
##                     ORDER[1],
##                     ANGLEUNIT["degree",0.0174532925199433]],
##                 AXIS["longitude",east,
##                     ORDER[2],
##                     ANGLEUNIT["degree",0.0174532925199433]],
##             ID["EPSG",4326]]],
##     ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
##         METHOD["Position Vector transformation (geog2D domain)",
##             ID["EPSG",9606]],
##         PARAMETER["X-axis translation",570.8,
##             ID["EPSG",8605]],
##         PARAMETER["Y-axis translation",85.7,
##             ID["EPSG",8606]],
##         PARAMETER["Z-axis translation",462.8,
##             ID["EPSG",8607]],
##         PARAMETER["X-axis rotation",4.998,
##             ID["EPSG",8608]],
##         PARAMETER["Y-axis rotation",1.587,
##             ID["EPSG",8609]],
##         PARAMETER["Z-axis rotation",5.261,
##             ID["EPSG",8610]],
##         PARAMETER["Scale difference",1.00000356,
##             ID["EPSG",8611]]]]
# transformace podle jiné existující vrstvy
brdy_krovak<-st_transform(brdy, st_crs(brdy_5514))
st_crs(brdy_krovak) # jaký crs má vrstva brdy_5514 ?
## Coordinate Reference System:
##   User input: +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 
##   wkt:
## BOUNDCRS[
##     SOURCECRS[
##         PROJCRS["unknown",
##             BASEGEOGCRS["unknown",
##                 DATUM["Unknown based on Bessel 1841 ellipsoid",
##                     ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
##                         LENGTHUNIT["metre",1,
##                             ID["EPSG",9001]]]],
##                 PRIMEM["Greenwich",0,
##                     ANGLEUNIT["degree",0.0174532925199433,
##                         ID["EPSG",9122]]]],
##             CONVERSION["unknown",
##                 METHOD["Krovak (North Orientated)",
##                     ID["EPSG",1041]],
##                 PARAMETER["Latitude of projection centre",49.5,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8811]],
##                 PARAMETER["Longitude of origin",24.8333333333333,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8833]],
##                 PARAMETER["Co-latitude of cone axis",30.2881397222222,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",1036]],
##                 PARAMETER["Latitude of pseudo standard parallel",78.5,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8818]],
##                 PARAMETER["Scale factor on pseudo standard parallel",0.9999,
##                     SCALEUNIT["unity",1],
##                     ID["EPSG",8819]],
##                 PARAMETER["False easting",0,
##                     LENGTHUNIT["metre",1],
##                     ID["EPSG",8806]],
##                 PARAMETER["False northing",0,
##                     LENGTHUNIT["metre",1],
##                     ID["EPSG",8807]]],
##             CS[Cartesian,2],
##                 AXIS["(E)",east,
##                     ORDER[1],
##                     LENGTHUNIT["metre",1,
##                         ID["EPSG",9001]]],
##                 AXIS["(N)",north,
##                     ORDER[2],
##                     LENGTHUNIT["metre",1,
##                         ID["EPSG",9001]]]]],
##     TARGETCRS[
##         GEOGCRS["WGS 84",
##             DATUM["World Geodetic System 1984",
##                 ELLIPSOID["WGS 84",6378137,298.257223563,
##                     LENGTHUNIT["metre",1]]],
##             PRIMEM["Greenwich",0,
##                 ANGLEUNIT["degree",0.0174532925199433]],
##             CS[ellipsoidal,2],
##                 AXIS["latitude",north,
##                     ORDER[1],
##                     ANGLEUNIT["degree",0.0174532925199433]],
##                 AXIS["longitude",east,
##                     ORDER[2],
##                     ANGLEUNIT["degree",0.0174532925199433]],
##             ID["EPSG",4326]]],
##     ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
##         METHOD["Position Vector transformation (geog2D domain)",
##             ID["EPSG",9606]],
##         PARAMETER["X-axis translation",570.8,
##             ID["EPSG",8605]],
##         PARAMETER["Y-axis translation",85.7,
##             ID["EPSG",8606]],
##         PARAMETER["Z-axis translation",462.8,
##             ID["EPSG",8607]],
##         PARAMETER["X-axis rotation",4.998,
##             ID["EPSG",8608]],
##         PARAMETER["Y-axis rotation",1.587,
##             ID["EPSG",8609]],
##         PARAMETER["Z-axis rotation",5.261,
##             ID["EPSG",8610]],
##         PARAMETER["Scale difference",1.00000356,
##             ID["EPSG",8611]]]]
# kde jsou ty brdy? no nevidím je protože to má jiný CRS
plot(st_geometry(hrcr),main="Kde jsou ty Brdy?") 
plot(st_geometry(brdy_32633),add=T) 

# musíme tedy transformovat jedno nebo druhé
hrcr_32633<-st_transform(hrcr, 32633)
plot(st_geometry(hrcr_32633),main="No jo, už máme správný CRS")
plot(st_geometry(brdy_32633),add=T,col="red")

Načítání rastrových GIS dat do R

# načítání jdnoduchou funkcí "raster"
DEM<-raster(paste0(cesta,"\\DEM_Jested_100m.tif"))
# kontrolí obrázek
# data tu jsou
plot(DEM)

# ale leží na správném místě? 
# ukáže crs vrstvy
DEM@crs
## CRS arguments:
##  +proj=krovak +axis=swu +lat_0=49.5 +lon_0=24.8333333333333 +alpha=0
## +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs
# zkušeným okem, nebo podle errorů při další alanýze  poznáte, 
# že nemá dobře definovaný crs, je  to křovák, ale šptně definovaný.
# tak mu řekneme, jaká je korektní definice křováka
# http://freegis.fsv.cvut.cz/gwiki/S-JTSK

DEM@crs<-crs("+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56")

# i rastr je možné vizualizovat interaktivně
# definujeme barvy
pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(DEM),
                    na.color = "transparent")

leaflet(DEM) %>% 
   addTiles() %>%  
   addRasterImage(DEM,colors = pal, opacity = 0.8) %>%
   addLegend(pal = pal, values = values(DEM),
             title = "Elevation [m]")

Příklad vektorové analýzy

# plocha polygonů
st_area(brdy)
## 345041249 [m^2]
# délka linií (obvod polygonů)
st_length(brdy)
## 195793.8 [m]
# buffer 5 km
brdy_5000<-st_buffer(brdy_32633,5000)
plot(st_geometry(brdy_5000), main="buffer 5 km")
plot(st_geometry(brdy_32633),add=T,col="red")

# centroidy
brdy_cen<-st_centroid(brdy_32633) 
## Warning in st_centroid.sf(brdy_32633): st_centroid assumes attributes are
## constant over geometries of x
plot(st_geometry(brdy_32633), main="centroid",graticule=T,axes=T)
plot(st_geometry(brdy_cen),add=T,col="blue",pch=19)

# průnik
chmu_brdy<-st_intersection(brdy,chmu_sf)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
plot(st_geometry(brdy), main="intersection - chmu stanice v brdech",graticule=T,axes=T)
plot(st_geometry(chmu_brdy),add=T,col="red")

Hromadné zpracování dat

Aneb co by nám v QGIS trvalo klikat hodiny můžeme v R naprogramovat během minut

# načteme všechny cesty k shp vrstvám v daném adresáři
parky.files<-list.files(path=paste0(cesta,"\\parky2"),pattern="*.shp$",full.names = T)

## když chceme všechny soubor do jednoho objektu
# načte první shpfile
parky.init<-st_read(parky.files[1],quiet = TRUE)

# připojí všechny další shp file do jendoho objektu  
for (i in 2:length(parky.files)) {
  parky.next<-st_read(parky.files[i],quiet = TRUE)
  parky.init<-rbind(parky.init,parky.next)
}  

# připojit metadat GIS join
meta.path<-paste0(cesta,"\\metadata_parky.csv")
metadata<-read.table(meta.path,sep=";",header=T,stringsAsFactors = T)
parky<-merge(parky.init,metadata)


# kontrolní obrázek
plot(st_geometry(hrcr))
plot(st_geometry(parky),add=T,col="green")

# interaktivně
distCol<- colorFactor(distinctColorPalette(32), parky.init$OBJECTID)
leaflet(parky) %>% 
  addTiles() %>%  
  addPolygons(color = ~distCol(OBJECTID),
              stroke = FALSE, fillOpacity = 0.8,
              label = paste0(parky$KAT,"_",parky$NAZEV))

Rastrové analýzy

Velice dobrá knihovna pro rastrové analýzy v R je whitebox, což je implementace jinak v Python/Rust napsaných geo algoritmů.

## instalace balíčku whitebox pro rastrové analýzy
# install.packages("whitebox", repos="http://R-Forge.R-project.org") # instalace
# whitebox::wbt_init() # inicializace
library(whitebox)
print(wbt_help())
##  [1] "WhiteboxTools Help"                                                                                        
##  [2] ""                                                                                                          
##  [3] "The following commands are recognized:"                                                                    
##  [4] "--cd, --wd       Changes the working directory; used in conjunction with --run flag."                      
##  [5] "-h, --help       Prints help information."                                                                 
##  [6] "-l, --license    Prints the whitebox-tools license."                                                       
##  [7] "--listtools      Lists all available tools. Keywords may also be used, --listtools slope."                 
##  [8] "-r, --run        Runs a tool; used in conjuction with --wd flag; -r=\"LidarInfo\"."                        
##  [9] "--toolbox        Prints the toolbox associated with a tool; --toolbox=Slope."                              
## [10] "--toolhelp       Prints the help associated with a tool; --toolhelp=\"LidarInfo\"."                        
## [11] "--toolparameters Prints the parameters (in json form) for a specific tool; --toolparameters=\"LidarInfo\"."
## [12] "-v               Verbose mode. Without this flag, tool outputs will not be printed."                       
## [13] "--viewcode       Opens the source code of a tool in a web browser; --viewcode=\"LidarInfo\"."              
## [14] "--version        Prints the version information."                                                          
## [15] ""                                                                                                          
## [16] "Example Usage:"                                                                                            
## [17] ">> .\\whitebox-tools.exe -r=lidar_info --cd=\"\\path\\to\\data\\\" -i=input.las --vlr --geokeys"           
## [18] ""
# pozor na diakritiku a mezery v cestě. Pro WBT nesmí být
# cesta k rastrovým datům
dem <-("c:\\Users\\manma\\Downloads\\DEM_Jested_100m.tif")

## stínovaný reliéf
# cesta k budoucímu výsledku
output<- ("c:\\Users\\manma\\Downloads\\Hillshade_100m.tif")

# spustí algoitmus pro výpočet stinovaného reliéfu 
wbt_hillshade(dem, output, azimuth = 315, altitude = 30, 
                    zfactor = 1,verbose_mode = FALSE)
## [1] "hillshade - Elapsed Time (excluding I/O): 0.4s"
# načtu výsledek jako rastr
hill<-raster(output)
# načtu původní dem jako rastr
elev<-raster(dem)
# vytisknu oba obrázky pro kontrolu
plot(hill, col=grey(0:100/100), legend=FALSE,main="stínovaný reliéf")
plot(elev, col=rainbow(25, alpha=0.35), add=TRUE)

## Sklon svahu
output<- ("c:\\Users\\manma\\Downloads\\Slope_100m.tif")
wbt_slope(dem,output,zfactor = 1, verbose_mode = FALSE)
## [1] "slope - Elapsed Time (excluding I/O): 0.3s"
slp<-raster(output)
plot(hill, col=grey(0:100/100), legend=FALSE,main="Sklon svahu")
plot(slp, col=heat.colors(25, alpha=0.2), add=TRUE)

## The terrain ruggedness index (TRI)
output<- ("c:\\Users\\manma\\Downloads\\TRI_100m.tif")
wbt_ruggedness_index(dem,output,zfactor = 1, verbose_mode = FALSE)
## [1] "ruggedness_index - Elapsed Time (excluding I/O): 0.4s"
tri<-raster(output)
plot(hill, col=grey(0:100/100), legend=FALSE,main="terrain ruggedness index")
plot(tri, col=cm.colors(25, alpha=0.3), add=TRUE)